ELISA
2022-08-15
Load in libraries
library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(formattable)
library(ggpubr)
library(pzfx)
library(DT)
library(cowplot)
source("scripts/utils/utils.R")Load data
# load data
d <- read_excel("processed_data/ELISA.xlsx")Observed data visualization
Supplementary Figure 3. Impact of mRNA-1273 Prime-boost Interval on S2-P-specific Serum Binding IgG Antibody Titers Through 24 Weeks Post-boost
#-------------------------------------------------------------------------------
# Supplementary Figure 1
#-------------------------------------------------------------------------------
# color panel
panel0 <-
c("#c85979",
"#6980ce",
"#4aac8d",
"#b460bd",
"#cd5d39",
"#b49541",
"#4ab0aa")
dosepanel <- c("#cb6939", "#8264cb")
LLOQ <- 25 # lower limit of quantificatin
# trend over time
raw1 <-
d %>%
filter(interval != "PBS") %>%
mutate(interval = factor(interval)) %>%
mutate(
.,
interval = forcats::fct_recode(
interval,
"8 weeks" = "8",
"6 weeks" = "6",
"4 weeks" = "4",
"3 weeks" = "3",
"2 weeks" = "2",
"1 week" = "1",
"prime only" = "0"
)
) %>%
mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
mutate(
day = forcats::fct_recode(
factor(day),
"Day 56 \n (Prime)" = "56",
"Day 64 \n (1wk post-boost)" = "64",
"Day 73 \n (2wk post-boost)" = "73",
"Day 87 \n (4wk post-boost)" = "87",
"Day 115 \n (8wk post-boost)" = "115",
"Day 143 \n (12wk post-boost)" = "143",
"Day 171 \n (16wk post-boost)" = "171",
"Day 199 \n (20wk post-boost)" = "199",
"Day 227 \n (24wk post-boost)" = "226"
)
) %>%
mutate(dose = forcats::fct_recode(factor(dose), "1ug" = "1", "10ug" = "10")) %>%
ggplot(aes(
x = factor(day),
y = value,
color = factor(interval)
)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(position = position_jitterdodge(jitter.width = .2)) +
geom_hline(aes(yintercept = log10(LLOQ), lty = "LLOQ")) +
scale_linetype_manual(values = 2) +
facet_wrap(~ dose, nrow = 3) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = " ",
y = expression(log[10] ~ S2P ~ IgG ~ titer),
color = "",
linetype = ""
) +
scale_color_manual(values = panel0) +
scale_y_continuous(breaks = 1:7) +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 16),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16)
)
# animal-level line plot
raw2 <-
d %>%
filter(interval != "PBS") %>%
mutate(
interval = factor(interval),
animal_id = paste("animal", animal_id),
animal_id = factor(
animal_id,
levels = c(
"animal 1",
"animal 2",
"animal 3",
"animal 4",
"animal 5",
"animal 6",
"animal 7",
"animal 8",
"animal 9",
"animal 10"
)
)
) %>%
mutate(
.,
interval = forcats::fct_recode(
interval,
"8 weeks" = "8",
"6 weeks" = "6",
"4 weeks" = "4",
"3 weeks" = "3",
"2 weeks" = "2",
"1 weeks" = "1",
"prime only" = "0"
)
) %>%
mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
mutate(dose = forcats::fct_recode(factor(dose), "1ug" = "1", "10ug" = "10")) %>%
ggplot(aes(
x = day,
y = value,
color = factor(animal_id)
)) +
geom_line() +
geom_jitter(position = position_jitterdodge(jitter.width = .2)) +
geom_hline(aes(yintercept = log10(LLOQ), lty = "LLOQ")) +
scale_linetype_manual(values = 2) +
facet_grid(dose ~ interval) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = " ",
y = expression(log[10] ~ S2P ~ IgG ~ titer),
color = "",
linetype = ""
) +
scale_color_manual(
values = c(
"#c25dba",
"#72b24a",
"#7561cf",
"#ce9b44",
"#807dc3",
"#7a7d37",
"#4faad9",
"#c95c3f",
"#4bac88",
"#c8577b"
)
) +
scale_y_continuous(breaks = 1:7) +
scale_x_continuous(breaks = c(56 , 87, 115, 143, 171, 199, 227)) +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 16),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16)
)
cowplot::plot_grid(plotlist = list( raw1, raw2 ), labels = c("A", "B"), nrow = 2, ncol = 1)
#ggsave(filename = "figures_tables/ELISA/plot_raw.png", width = 16, height=16)Model fitting
Residual diagnostics when model day as a continuous variable
#-------------------------------------------------------------------------------
# Fit the model
#-------------------------------------------------------------------------------
dtemp <- filter(d, interval != "PBS" & interval != "0") %>%
mutate(
dose = factor(dose),
day = as.numeric(day),
animal_id = factor(animal_id),
interval = factor(interval)
) %>%
mutate(
animal_unique_id = paste0(interval, "_", animal_id),
animal_unique_id = factor(animal_unique_id)
)
# fit a gam model
gam_fit <-
mgcv::gam(
value ~ interval + dose + s(day, k = 9) + interval:day +
interval:dose + dose:day + dose:day:interval +
s(interval, bs = "re") + s(interval, day, bs = "re") ,
data = dtemp,
method = 'REML'
)
# Residual diagnostics figure
p <- assumption_check(gam_fit)
title <-
ggdraw() + draw_label(
"ELISA residual diagnostics when model day as a continuous variable",
x = 0.03,
hjust = 0
)
plot_grid(
title,
p,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
# fit a gam linear model
gam_fit_linear <-
mgcv::gam(
value ~ interval + dose + day + interval:day +
interval:dose + dose:day + dose:day:interval +
s(interval, bs = "re") + s(interval, day, bs = "re") ,
data = dtemp,
method = 'REML'
)
# fit a gam model at df = 8
gam_fit_df8 <-
mgcv::gam(
value ~ interval + dose + s(day, k = 8) + interval:day +
interval:dose + dose:day + dose:day:interval +
s(interval, bs = "re") + s(interval, day, bs = "re") ,
data = dtemp,
method = 'REML'
)Diagnostics
The model without smoothing spline on day does not fit
as good as the selected model, where the selected model has AIC = 883
and BIC = 1042 and the model without smoothing spline has AIC = 1640 and
BIC = 1764.
The selected model has degrees of freedom equals 9 for the smoothing spline, which is the maximum degrees of freedom allowed given the unique covariate combinations. AIC and BIC tests suggested degree of freedom equals 9 has the best fit. The model has df = 8 smoothing spline has AIC = 984 and BIC = 1138.
Therefore, we add smoothing spline with df = 9 to variable
day.
# Residual diagnostics figure
p <- assumption_check(gam_fit_linear)
title <-
ggdraw() + draw_label(
"ELISA residual diagnostics when model day as a continuous variable",
x = 0.03,
hjust = 0
)
plot_grid(
title,
p,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
Figure 2. S2-P-specific Serum Binding IgG Antibody Titers
#-------------------------------------------------------------------------------
# Figure 2. Spike-specific Serum Binding IgG Antibody Titers
#-------------------------------------------------------------------------------
plot1 <- ggemmeans(gam_fit, terms = c("day", "interval", "dose"))
plot(plot1) +
scale_color_manual(values = panel0) +
scale_x_continuous(
breaks = c(56 , 64, 73, 87, 115, 143, 171, 199, 226),
labels = c(
"56(prime)",
"64(1wk post-boost)",
"73(2wk post-boost)",
"87(4wk post-boost)",
"115(8wk post-boost)",
"143(12wk post-boost)",
"171(16wk post-boost)",
"199(20wk post-boost)",
"227(24wk post-boost)"
)
) +
labs(
title = expression(Figure2A ~ Predicted ~ values ~ "for" ~ log[10] ~ S2P ~
IgG ~ titer),
y = expression(log[10] ~ S2P ~ IgG ~ titer)
) +
theme(panel.grid.major = element_line(colour = "grey100")) +
theme(
axis.text = element_text(size = 12),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
axis.title = element_text(size = 18),
strip.text = element_text(size = 18),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
plot.title = element_text(size = 20),
legend.position = "bottom"
) +
facet_grid( ~ facet) 
# plot effect of dose treating dose as numerical variable
plot2 <-
ggemmeans(
gam_fit,
terms = c("day", "dose", "interval"),
condition = c("day", "interval")
)
plot(plot2) +
#scale_color_manual(values=dosepanel) +
scale_x_continuous(
breaks = c(56 , 64, 73, 87, 115, 143, 171, 199, 226),
labels = c(
"56(prime)",
"64(1wk post-boost)",
"73(2wk post-boost)",
"87(4wk post-boost)",
"115(8wk post-boost)",
"143(12wk post-boost)",
"171(16wk post-boost)",
"199(20wk post-boost)",
"226(24wk post-boost)"
)
) +
labs(
title = expression(Figure2B ~ Predicted ~ values ~ "for" ~ log[10] ~ S2P ~
IgG ~ titer),
y = expression(log[10] ~ S2P ~ IgG ~ titer)
) +
theme(panel.grid.major = element_line(colour = "grey100")) +
theme(axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
legend.position = "bottom") +
facet_wrap( ~ facet)
Residual diagnostics when model day as a factor
# fit the model for mean comparison
dtemp <- filter(d, interval != "PBS" & interval != "0") %>%
mutate(
dose = as.factor(dose),
animal_id = factor(animal_id),
interval = factor(interval),
day = factor(day)
) %>%
mutate(animal_id = paste0(interval, "_", animal_id)) %>%
mutate(animal_id = paste0(animal_id , "_", dose))
lm_fit <-
lmer(value ~ interval * dose * day + (1 | animal_id) , data = dtemp)
p2 <- assumption_check(lm_fit)
# Residual diagnostics figure
title <-
ggdraw() + draw_label(
"ELISA residual diagnostics when model day as a discrete variable",
x = 0.03,
hjust = 0
)
plot_grid(
title,
p2,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
Supplementary Figure 1. Predicted Fold-change of Spike-specific Serum Binding IgG Antibody Titers
# The results indicate that the peak immune response/ titer levels are achieved at 2 weeks post prime
grid <- ref_grid(lm_fit, cov.keep = c("day", "dose", "interval"))
emm_fit_day <- emmeans(grid, ~ day | interval + dose)
#-------------------------------------------------------------------------------
# Supplementary Figure 2. Predicted Fold-change of Spike-specific Serum Binding IgG Antibody Titers.
#-------------------------------------------------------------------------------
set.seed(101)
summary(
contrast(
emm_fit_day,
list(
"Day 73 vs Day 56 \n(2wk post-boost) vs (Prime)" = c(-1, 0, 1, 0, 0, 0, 0, 0, 0),
"Day 227 vs Day 56 \n(24wk post-boost) vs (Prime)" = c(-1, 0, 0, 0, 0, 0, 0, 0, 1)
)
),
infer = TRUE,
adjust = "mvt",
level = .95
) %>%
as_tibble() %>%
mutate(
.,
interval = forcats::fct_recode(
interval,
"8 weeks" = "8",
"6 weeks" = "6",
"4 weeks" = "4",
"3 weeks" = "3",
"2 weeks" = "2",
"1 weeks" = "1"
)
) %>%
mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
mutate(dose = forcats::fct_recode(factor(dose), "1ug" = "1", "10ug" = "10")) %>%
ggplot(aes(x = contrast , y = estimate, colour = interval)) +
geom_point(position = position_dodge(width = 0.6), size = 2) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 0.6),
size = 1
) +
scale_color_manual(values = panel0) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.5, 1, 3, 10, 30, 100)), labels = c(0.5, 1, 3, 10, 30, 100)) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
labs(x = "",
y = expression(Estimated ~ fold ~ change ~ of ~ S2P ~ IgG ~ titer)) +
facet_grid( ~ dose) +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 18),
strip.text = element_text(size = 18),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
plot.title = element_text(size = 20),
legend.position = "bottom"
)
Supplementary Table 1. Statistical comparisons of S2-P IgG antibody titer in mice between different mRNA-1273 dosing intervals.
#-------------------------------------------------------------------------------
# Supplementary Table 1
#-------------------------------------------------------------------------------
grid <- ref_grid(gam_fit, cov.keep = c("day", "dose", "interval"))
em_fit <- emmeans::emmeans(grid, ~ interval | day + dose)
set.seed(100)
sum_fit <-
summary(
pairs(em_fit , reverse = TRUE),
infer = TRUE,
adjust = "mvt",
level = .95
)
# generate tables
sum_fit %>%
as_tibble() %>%
dplyr::select(contrast, estimate, day, dose, p.value) %>%
mutate(estimate = 10 ^ estimate) %>%
mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value"
) %>%
filter(`Adjusted P-value` < 0.05) %>%
DT::datatable(caption = "Statistical comparisons of S2-P IgG antibody titer in mice between different mRNA-1273 dosing intervals") %>%
formatSignif(columns = c('Fold change', 'Adjusted P-value'),
digits = 7)Supplementary Figure 4. Statistical comparisons of S2-P IgG Antibody Titers between different mRNA-1273 dosing intervals(Supplementary Table 1 visualization)
sum_fit %>%
as_tibble() %>%
mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
mutate(day = as.factor(day)) %>%
mutate(
day = forcats::fct_recode(
factor(day),
"Day 56 \n(Prime)" = "56",
"Day 64 \n(1wk post-boost)" = "64",
"Day 73 \n(2wk post-boost)" = "73",
"Day 87 \n(4wk post-boost)" = "87",
"Day 115 \n(8wk post-boost)" = "115",
"Day 143 \n(12wk post-boost)" = "143",
"Day 171 \n(16wk post-boost)" = "171",
"Day 199 \n(20wk post-boost)" = "199",
"Day 227 \n(24wk post-boost)" = "226"
)
) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
"Day" = "day"
) %>%
mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
mutate(dose = forcats::fct_recode(factor(dose), "1ug" = "1", "10ug" = "10")) %>%
ggplot(aes(x = `Pairwise comparison` , y = `Fold change`, color = Day)) +
geom_point(position = position_dodge(width = 2), size = 3) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 2),
size = 1
) +
scale_color_manual(
values = c(
"#ABCFE5",
"#93C4DE",
"#78B5D8",
"#60A6D1",
"#4A97C9",
"#3787C0",
"#2575B7",
"#1664AB",
"#09539D",
"#084185"
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.5, 1, 3, 10, 30, 100)),
labels = c(0.5, 1, 3, 10, 30, 100)) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
labs(x = "",
y = expression(Estimated ~ fold ~ change ~ of ~ S2P ~ IgG ~ titer)) +
facet_grid(`Pairwise comparison` ~ dose,
scales = "free",
space = "free") +
coord_flip() +
theme(
axis.text = element_text(size = 18),
axis.title = element_text(size = 18),
strip.text = element_text(size = 18),
legend.text = element_text(size = 16,
margin = margin(t = 10)),
legend.title = element_text(size = 18),
plot.title = element_text(size = 20),
strip.text.y = element_blank(),
legend.position = "right"
)
ggsave(filename = "figures_tables/ELISA/ELISA_emmeans.png", width = 16, height=18)Session information
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.1.1 DT_0.24 pzfx_0.3.0 ggpubr_0.4.0
## [5] formattable_0.2.1 lme4_1.1-28 Matrix_1.2-18 ggeffects_1.1.3
## [9] emmeans_1.7.5 mgcv_1.8-31 nlme_3.1-144 readxl_1.4.0
## [13] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4
## [17] readr_2.1.2 tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6
## [21] tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] pbkrtest_0.5.1 fs_1.5.2 lubridate_1.8.0
## [4] insight_0.18.2 httr_1.4.3 tools_3.6.3
## [7] backports_1.4.1 bslib_0.4.0 sjlabelled_1.2.0
## [10] utf8_1.2.2 R6_2.5.1 DBI_1.1.3
## [13] colorspace_2.0-3 withr_2.5.0 tidyselect_1.1.2
## [16] compiler_3.6.3 textshaping_0.3.6 cli_3.3.0
## [19] rvest_1.0.2 xml2_1.3.3 labeling_0.4.2
## [22] sass_0.4.2 scales_1.2.0 mvtnorm_1.1-3
## [25] systemfonts_1.0.2 digest_0.6.29 minqa_1.2.4
## [28] rmarkdown_2.14 pkgconfig_2.0.3 htmltools_0.5.3
## [31] highr_0.9 dbplyr_2.2.1 fastmap_1.1.0
## [34] htmlwidgets_1.5.4 rlang_1.0.4 rstudioapi_0.13
## [37] farver_2.1.1 jquerylib_0.1.4 generics_0.1.3
## [40] jsonlite_1.8.0 crosstalk_1.2.0 car_3.1-0
## [43] googlesheets4_1.0.1 magrittr_2.0.3 Rcpp_1.0.9
## [46] munsell_0.5.0 fansi_1.0.3 abind_1.4-5
## [49] lifecycle_1.0.1 stringi_1.7.8 yaml_2.3.5
## [52] carData_3.0-5 MASS_7.3-51.5 grid_3.6.3
## [55] parallel_3.6.3 crayon_1.5.1 lattice_0.20-38
## [58] haven_2.5.0 splines_3.6.3 hms_1.1.1
## [61] knitr_1.39 pillar_1.8.0 boot_1.3-24
## [64] estimability_1.4.1 ggsignif_0.6.3 reprex_2.0.1
## [67] glue_1.6.2 evaluate_0.16 modelr_0.1.8
## [70] vctrs_0.4.1 nloptr_1.2.2.3 tzdb_0.3.0
## [73] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [76] cachem_1.0.6 xfun_0.32 xtable_1.8-4
## [79] broom_1.0.0 coda_0.19-4 rstatix_0.7.0
## [82] ragg_1.1.3 googledrive_2.0.0 gargle_1.2.0
## [85] ellipsis_0.3.2